Estudiante Ing. Forestal
NASA
“Desde 1880, el nivel del mar global ha aumentado 20 centimetros. Para el 2100 se proyecta que aumente entre 30 y 122 centimetros más.”
Fuentes
pacman::p_load(TSstudio, tsbox, tsoutliers, tidyverse, lubridate, xts, magrittr, forecast, dygraphs) url <- 'https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_mm_mlo.txt' cov <- read.table(file = url) %>% select(3, 4) %>% rename(fecha = 1, co2 = 2) %>% mutate(fecha = date_decimal(fecha) %>% floor_date(unit = 'month')) %>% filter(fecha <= '2013-12-01') data <- read.csv('gmsl.csv') %>% mutate(fecha = date_decimal(Time) %>% floor_date(unit = 'month')) %>% rename(gsml = 2) %>% select(4, 2)%>% filter(fecha > '1958-02-01') %>% left_join(y = cov) %>% {xts(x = .[,2:3], order.by = .$fecha)} rm(cov) dygraph(data, main = 'Nivel medio del mar global vs Concentración de CO2 en la atmósfera') %>% dySeries("co2", axis = 'y2', label = 'CO2') %>% dySeries('gsml', label = 'GSML') %>% dyAxis('y', label = 'Nivel del mar (mm)') %>% dyAxis('y2', label = 'CO2 (ppm)', independentTicks = T) %>% dyRangeSelector()
ts_split(ts.obj = ts_ts(data), sample.out = 12) %T>% {assign(x = "train", value = .$train, envir = .GlobalEnv)} %T>% {assign(x = "test", value = .$test, envir = .GlobalEnv)}
acf(train[,'gsml'], main = '')
tseries::adf.test(x = train[,'gsml'])
Augmented Dickey-Fuller Test data: train[, "gsml"] Dickey-Fuller = -2.7778, Lag order = 8, p-value = 0.2491 alternative hypothesis: stationary
matriz_diseño <- model.matrix(~ -1 + train[,'co2']) colnames(matriz_diseño) <- 'co2' modelo <- Arima(y = train[,'gsml'], order = c(2,1,3), xreg = matriz_diseño, include.drift = T) modelo %>% lmtest::coeftest()
z test of coefficients: Estimate Std. Error z value Pr(>|z|) ar1 0.797070 0.054991 14.4947 < 2.2e-16 *** ar2 -0.095533 0.054919 -1.7395 0.08194 . ma1 -0.339090 0.040237 -8.4274 < 2.2e-16 *** ma2 0.208116 0.038717 5.3753 7.646e-08 *** ma3 -0.693544 0.028246 -24.5535 < 2.2e-16 *** drift 0.208518 0.041471 5.0281 4.955e-07 *** co2 -0.171228 0.087716 -1.9521 0.05093 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
equatiomatic::extract_eq(modelo, wrap = T, use_coefs = T)
\[ \begin{alignat}{2} &\text{let}\quad &&y_{t} = \operatorname{y}_{\operatorname{0}} +0.21\operatorname{t} -0.17\operatorname{co2}_{\operatorname{t}} +\eta_{t} \\ &\text{where}\quad &&(1 -0.8\operatorname{B} +0.1\operatorname{B}^{\operatorname{2}} )\ (1 - \operatorname{B}) \eta_{t} \\ & &&= (1 -0.34\operatorname{B} +0.21\operatorname{B}^{\operatorname{2}} -0.69\operatorname{B}^{\operatorname{3}} )\ \varepsilon_{t} \\ &\text{where}\quad &&\varepsilon_{t} \sim{WN(0, \sigma^{2})} \end{alignat} \]
modelo %>% checkresiduals(theme = theme_bw(), test = F)
modelo %>% residuals() %>% shapiro.test()
Shapiro-Wilk normality test data: . W = 0.99861, p-value = 0.8921
modelo %>% residuals() %>% Box.test(type = 'Ljung-Box')
Box-Ljung test data: . X-squared = 0.029686, df = 1, p-value = 0.8632
tso(y = train[,'gsml'], xreg = train[,'co2'])
Series: Regression with ARIMA(1,1,3) errors Coefficients: ar1 ma1 ma2 ma3 xreg 0.7152 -0.2650 0.1905 -0.6741 -0.1222 s.e. 0.0439 0.0387 0.0333 0.0275 0.0851 sigma^2 = 3.074: log likelihood = -1299.53 AIC=2611.06 AICc=2611.19 BIC=2637.99 No outliers were detected.
prediccion <- forecast(object = modelo, h = 12, level = 95, xreg = test[,'co2']) serie <- prediccion %>% as_tibble() %>% ts(start = c(2013, 1), end = c(2013, 12), frequency = 12) %>% {cbind(data[, 'gsml'], .)} colnames(serie) <- c('gsml', 'predicho', 'lwr', 'upr') serie %>% dygraph() %>% dySeries('gsml', label = 'Observada') %>% dySeries(c('lwr', 'predicho', 'upr'), label = 'Predicho') %>% dyRangeSelector()
accuracy(prediccion, test[,'gsml'])[,c(2,5)]
RMSE MAPE Training set 1.719463 Inf Test set 9.180928 12.04332
